

To model evolving material configurations with two displacment fields, we will specify the following through defining user functions\-: \par

\begin{DoxyItemize}
\item Boundary conditions \par

\item Load stepping \par

\item Derived fields for output \par

\item Constitutive model \par

\item Parameter values \par

\item Weak form of the P\-D\-E \par

\end{DoxyItemize}

See the paper \char`\"{}\-A variational treatment of material configurations with application to interface motion and microstructural evolution\char`\"{}, G. Teichert, et al. (Journal of the Mechanics and Physics of Solids, 2017).

First, we include the header file declaring the required user functions. These functions will be defined in this file.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Now, we first define any optional user functions. Optional user functions have a default definition that can be redefined by the user using a function pointer. This will be done in the {\ttfamily define\-Parameters} function. The available list of optional user functions includes\-: {\ttfamily boundary\-Conditions}, {\ttfamily scalar\-Initial\-Conditions}, {\ttfamily vector\-Initial\-Conditions}, {\ttfamily load\-Step}, {\ttfamily adaptive\-Time\-Step}, and {\ttfamily project\-Fields}. In this example, we redefine the {\ttfamily boundary\-Conditions} and {\ttfamily project\-Fields} functions, while using the default functions for the others.

{\bfseries  The {\ttfamily boundary\-Conditions} function }

This function defines Dirichlet boundary conditions using Pet\-I\-G\-A's {\ttfamily I\-G\-A\-Set\-Boundary\-Value} function. The arguments to this function are as follows\-: the iga object (user.\-iga), the \char`\"{}axis\char`\"{} (0, 1, or 2, corresponding to the x, y, or z-\/axis), the \char`\"{}side\char`\"{} (0 or 1), the \char`\"{}dof\char`\"{}, and the \char`\"{}value\char`\"{} that is to be imposed. Note that this can only set a uniform value for a degree-\/of-\/freedom on any side. Here, we fix all displacements for both configurational and total displacement fields at x=0 (axis=0,side=0). We also define vertical configurational and total displacements at x=10 (axis=0,side=1). Note that the magnitude of the displacement is dependent on the {\ttfamily scale} parameter, which is defined in the next function, {\ttfamily load\-Step}.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily load\-Step} function }

This function allows us to update the {\ttfamily scale} parameter and call the {\ttfamily boundary\-Conditions} function with the updated value at every time step. This is useful when the Dirichlet boundary condition is too large to allow convergence when applied all at once. Here, the value for /scale is taken from the current \char`\"{}time\char`\"{}.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily project\-Fields} function }

If there are field values derived from the solution fields that are of interest, we can compute these values at each quadrature point and project the value to the nodes. Here, we compute the value {\ttfamily eta2}, which is a function of the configurational displacement that reflects the local rectangular variant of the crystal structure.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily constitutive\-Model} function }

This function defines the standard and configurational stresses, strain energy densities, and kinematic values that appear in the resiudal. Note that it is used only in this file (by the {\ttfamily residual} functions), so it is not a class member function nor does it have an associated function pointer.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


The first part of this function computes variables associated with the standard strain energy density function $\psi^\mathrm{S}$\-:

\begin{eqnarray*} \boldsymbol{E} &=& \frac{1}{2}(\boldsymbol{F}^T\boldsymbol{F} - \mathbbm{1})\\ \alpha_I(\boldsymbol{\chi}) &=& \alpha\Lambda_I(\boldsymbol{\chi}), \, \Lambda_I^2 = \sum_{i=1}^3\chi_{iI}^2\\ \mathbb{C}(\boldsymbol{\chi}) &=& \beta(\mathbbm{1}\otimes\mathbbm{1}) + 2\mu\mathbb{I} + \sum_{I=1}^3 \left(\alpha_I(\boldsymbol{\chi}) -\beta - 2\mu \right) \boldsymbol{e}_I\otimes \boldsymbol{e}_I\otimes\boldsymbol{e}_I\otimes \boldsymbol{e}_I\\ \psi^\mathrm{S} &=& \frac{1}{2} \boldsymbol{E}:\mathbb{C}(\boldsymbol{\chi}):\boldsymbol{E} \\ &=& \frac{1}{2}\left[\beta\mathrm{tr}(\boldsymbol{E})^2 + 2\mu(\boldsymbol{E}:\boldsymbol{E}) + \sum_{I=1}^3 E_{II}^2\left(\alpha_I(\boldsymbol{\chi}) -\beta - 2\mu \right) \right]\\ \frac{\partial \psi^\mathrm{S}}{\partial \boldsymbol{\chi}} &=& \frac{1}{2} \sum_{I=1}^3 \frac{\alpha}{\Lambda_I}\boldsymbol{\chi}(\boldsymbol{e}_I\otimes \boldsymbol{e}_I) E_{II}^2 \end{eqnarray*}


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


The second part computes variables associated with the material strain energy density function $\psi^\mathrm{M}$\-:

\begin{eqnarray*} \boldsymbol{\Theta} &=& \half(\boldsymbol{\chi}^T\boldsymbol{\chi} - \mathbbm{1})\\ \eta_1 &=& \Theta_{11} + \Theta_{22},\, \eta_2 = \Theta_{11} - \Theta_{22},\, \eta_6 = \Theta_{12}\\ \psi^\mathrm{M} &=& \frac{d}{s^2}\left(\eta_1^2 + \eta_6^2\right) -\frac{2d}{s^2}\eta_2^2 + \frac{d}{s^4}\eta_2^4 + \frac{l^2d}{s^2}|\nabla^0 \eta_2|^2\\ \boldsymbol{D} &=& \frac{\partial \psi^\mathrm{M}}{\partial \boldsymbol{\chi}},\, \boldsymbol{B} = \frac{\partial \psi^\mathrm{M}}{\partial \nabla^0\boldsymbol{\chi}}\\ \end{eqnarray*}


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily define\-Parameters} function }

The user is required to define the {\ttfamily define\-Parameters} and {\ttfamily residual} functions. The {\ttfamily define\-Parameters} defines variables and functions in the {\ttfamily App\-Ctx} object. The {\ttfamily App\-Ctx} object is defined in the app\-Ctx.\-h file. This function is used to define any values in {\ttfamily user} that will be needed in the problem. It is also used to set any function pointers for user functions that we have redefined.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We define the mesh, domain, and total applied displacement.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We specify the number of vector and scalar solution and projection fields by adding the name of each field to their respective vector. Here, we have two vector solution fields (the configurational displacement and total displacement) and one scalar projection field (eta2). We do not use any scalar solution or vector projection fields in this example.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We can specify the polynomial order of the basis splines, as well as the global continuity. Note that the global continuity must be less than the polynomial order. Here, we use a quadratic basis function with C-\/1 continuity because of the gradient elasticity terms.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


We specify the value {\ttfamily dt\-Val}, which is a load step in this problem, as well as the {\ttfamily total\-Time}, or total loading. We can also specify a location (iteration number and time) to restart.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we redirect the desired user function pointers to the {\ttfamily boundary\-Conditions}, {\ttfamily project\-Fields}, and {\ttfamily load\-Step} functions that we defined above. This completes the {\ttfamily define\-Parameters} function.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


{\bfseries  The {\ttfamily residual} function }

The residual function defines the residual that is to be driven to zero. This is the central function of the code. It is set up to follow the analytical weak form of the P\-D\-E. It has a number of arguments that give problem information at the current quadrature point (see Example 1 and the documentation for the arguments' classes for further information).

The example code here implements the following weak form (neglecting body force and traction boundary conditions)\-:

\begin{eqnarray*} 0 &=& \int \limits_{{\Omega}_0} \nabla^0\nabla^0\boldsymbol{W}\,\vdots\,\boldsymbol{B} \,\mathrm{d}V_0+\\ &\phantom{=}& +\int \limits_{{\Omega}_0} \nabla^0\boldsymbol{W}:\left[\boldsymbol{D}+J_\chi\left(\mathcal{E} \boldsymbol{\chi}^{-T} + \frac{\partial \psi^\mathrm{S}}{\partial \boldsymbol{\chi}} \right) \right] \, \mathrm{d}V_0\\ &\phantom{=}& + \int \limits_{{\Omega}_0} \nabla^0\bar{\boldsybmol{w}}:\left(J_\chi\boldsymbol{P} \boldsymbol{\chi}^{-T} \right) \, \mathrm{d}V_0 \end{eqnarray*}

After calling the {\ttfamily constitutive\-Model} function to compute the necessary values, we compute the residual in a manner very similar to the analytical weak form\-:


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


Finally, we include a file that instatiates the template functions {\ttfamily define\-Parameters} and {\ttfamily residual}. This bit of code will generally be the same for any problem (unless you decide to use a different automatic differentation library), the user does not need to modify it.


\begin{DoxyCodeInclude}

\end{DoxyCodeInclude}


\section*{The complete code }


\begin{DoxyCodeInclude}
\end{DoxyCodeInclude}
 